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TECHNICAL  NOTE  1965-21 

14  JUNE  1965 


LEXINGTON 


MASSACHUSETTS 


A  COMPARISON  OF  SEISMIC  ARRAY  PROCESSING  SCHEMES 


E.  J.  Kelly 


ABSTRACT 

It  is  our  purpose  in  this  note  to  discuss  three  of  the  many  approaches  to  seismic 
array  processing  from  the  theoretical  point  of  view.  The  three  are:  1)  maximum  - 
likelihood  processing,  2)  the  minimum-variance,  unbiased  estimator  (MVU)  approach 
used  by  Levin,  and  3)  multichannel  Wiener  filtering.  A  feature  common  to  these 
techniques  is  the  formation  of  a  single  output  waveform  which  serves  as  an  estimator 
of  the  unknown  signal  coming  from  a  fixed  direction.  We  refer  to  such  an  output  as  a 
"beam.  ” 

It  will  be  shown  that  1)  reduces  to  2)  in  the  case  of  gaussian  noise  and  known 
signal  parameters  (i.  e.  ,  known  epicenter)  and  that  3)  is  related  very  simply  to  2). 

In  fact,  the  Wiener  filtering  output  can  be  obtained  directly  from  the  MVU  output.  We 
treat  both  the  sampled-data  and  continuous -time  cases.  The  note  contains  no  new 
results. 
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I.  Introduction 


Suppose  that  the  array  seismometers  are  identical  (say  short -period  vertical 
instruments  and  are  situated  in  a  horizontal  plane  at  the  points  r  k,  k  =  1, . . .  ,N.  We 
write  the  output,  x^(t),  of  the  k*-^1  seismometer  as  a  sum  of  signal  and  noise: 

x^t)  =  sk(t)  +  nk(t)  (1) 

We  suppose  that  the  noise  components  have  zero  mean  and  are  (wide -sense)  stationary 
with  covariance  matrix 

=  E  ^(t)  n^(t’)  (2) 

(Here  E  stands  for  expectation  or  ensemble  average. )  The  signal  is  an  unknown  wave¬ 
form,  propagating  across  the  array  with  horizontal  phase  velocity  v  and  coming  from 
a  direction  bearing  0  to  the  East  of  North.  If  this  signal  has  the  form  S(t)  at  the  origin 
of  coordinates,  it  has  the  form  S(t-a*r)  at  point  r,  where 

sin  0 

a  = 
x  v 

cos  0 

a  =  — — — 

y  v 

in  a  coordinate  system  where  the  x-axis  points  East  and  the  y-axis  points  North.  Thus, 

sk(t)  =  S(t-a*r  k)  (3) 

We  assume,  at  first,  that  the  signal  parameters  (ctx  and  a  ,  or  v  and  0)  are 
known,  which  is  equivalent  to  knowing  the  epicenter  if  the  depth  of  the  event  is  assumed 
or  known.  Then  the  array  can  be  "steered"  by  delaying  the  output  of  the  k^1  seismom¬ 
eter  by  the  amount  a-r^.  We  call  this  output  x^(t): 

x^(t)  =  Xj^t+a  •  r  R)  =  S(t)  +  n^(t)  (4) 

where 

n],(t)  =  n^t+a-TjP 
now  has  the  covariance  matrix 

R^(t-t’)  =  En^(t)nJ(t')  =  Rk/.[t-t’+a-(rk-r^)]  (5) 

The  x^(t)  are,  of  course,  equivalent  to  the  xk(t)  and  no  information  has  been  destroyed. 
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II.  The  Minimum -Variance  Unbiased  Estimator 


In  this  completely  ad  hoc  approach  we  postulate  a  linear  combination  of  the 
steered  array  outputs  which  is  to  provide  a  minimum  variance  unbiased  estimate  of 
the  unknown  signal  S(t).  Such  an  output  beam  has  the  form 

N 

y'(t)  =  Tj  J  h.(T)  X*  (t-T)  dT  (6) 

k=l  *  k 

i.  e. ,  the  steered  array  outputs  are  filtered  and  then  summed.  We  do  not  impose 
realizability  on  the  filter  functions,  hj,(t)f  and  all  time  integrals  run  from  — 00  to  +°°. 
Since 

E  x^(t)  =  S(t), 

the  expected  output  is 

N 

E  yfr)  =  I  J  Mt)  s  (t-T)  dT. 

k=l 

If  y*(t)  is  to  be  an  unbiased  estimator  of  S(t)  for  any  S(t),  we  must  have 

N 

Z  h  (t)  =  6(t)  .  (7) 

k=l  k 

The  variance  of  y’(t),  which  is  a  constant  since  y’(t)  is  a  stationary  process,  is 

o2  =  E  [y’(t)  -  S  (t)  ] 2 

N  2 

=  E  {  Z  I  Mt)  n!  (t-T)  dT  I 
k=l  k 

N 

=  Z  SS  hk(T)h^(T’)  R’  (T’-T)dTdT’  (8) 

k=l 


The  functions  h^(t)  are  to  be  chosen  to  minimize  expression  (8)  subject  to  constraint  (7). 
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In  the  sampled-data  version  of  the  problem,  we  suppose  that  the  steered  array 
outputs  are  sampled  every  6  seconds,  providing  a  sequence 


rk(i)  =  x^(is) 


i  =  0,  ±1,  ±2,  etc. 


and  the  output  beam  is  the  sequence 

N  v 

ri(i)  =  Z  Z  ek0)  ?k(i-j)  (9) 

k=l j=-v 

where  it  is  supposed  that  the  filters,  0^(i),  are  of  length  (2v+l)  samples,  symmetrically 
placed  about  zero  (this  last  is  non-essential  assumption).  Since 


E  ^(i)  =  £(i)  =  S(iS), 


the  analogue  of  (7)  is 


N 

Z 


k=l 


(10) 


In  terms  of  the  covariance  matrix 


■  vi6-j6)=  pko-‘>' 

the  output  variance  is  just 

N  v 

k,  1=1  i. 

We  solve  the  discrete  version  first,  using  the  calculus  of  variations  to  minimize 
(11)  subject  to  the  set  of  (2v+l)  constraints  (10).  Using  the  symmetry  of  the  covariance 
matrix,  we  find 


Z  P^O-i)  ek(i)  \d)  (11) 

j=-v 


o2  =  E  r\(if  =  Z 


N 


v 


for  the  variation  of  o3  under  variations,  6 ©^(i),  in  the  filter  weights.  We  multiply  the 
variation  of  the  im  constraint  by  the  Lagrange  multiplier  2  \(i),  sum  over  i  and  add 
the  result  to  6  c3.  The  6  0^(i)  are  now  treated  as  independent  variations,  which  leads 
immediately  to  the  set  of  linear  equations 

N  v 

2  2  0^(j)  +  Mi)  =  0  (12) 

•^=1  j=-v 

where  k  =  1, . . .  ,N  and  i  =  -v, . . .  ,+v. 

Equations  (12)  could  be  solved  for  0^(i)  in  terms  of  the  undetermined  k(i),  and 
the  latter  determined  by  substitution  in  the  equations  of  constraint.  However,  since  the 
constraints  are  also  linear  in  the  9^(i) »  we  may  adjoin  equations  (10)  to  equations  (12), 
treat  the  k(i)  as  unknowns,  and  solve  at  once  for  the  0^(i)  and  the  A.(i).  We  observe  that 
(12)  implies  that  o3  is  equal  to 

N  v 

o2  =  2  2  ek(i)[-Mi)] 

k=l  i=-v 

v 

=  -  2  6,  Mi)  =  -  Mo)  ,  (13) 

i=-v  1,0 

hence  the  value  of  k(o),  at  least,  is  of  interest  since  it  gives  us  the  actual  minimum 
noise  variance  attained.  If  we  define 

eN+1(i)  =  MO  (14) 


we  find  that  the  system  of  equations  can  be  written 


N+l  v 

2  .2  pfc£p_1)  e^i)  =  6k,N+i 

L=1  j=-v 


6. 

J,o 


where  we  have  extended  the  definition  of  the  matrix  P  as  follows: 


(15) 
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and 


The  set  of  (N+l)  (2v+l)  equations  (15)  is  equivalent*  to  the  set  given  by  Levin  in 
Reference  (1). 

In  the  continuous -time  case,  we  could  carry  out  a  similar  analysis  in  time  do¬ 
main,  by  introducing  a  Lagrangian  function,  \(t),  to  incorporate  the  infinity  of  constraints 
contained  in  (7),  and  leading  to  a  system  of  integral  equations  for  the  functions  hk(t)  and 
Mt)-  However,  if  neither  realizability  nor  finite  memory  is  imposed  on  the  hk(t),  the 
solution  is  more  easily  obtained  in  frequency  domain.  We  introduce  the  transform  filter 
functions 


Hk(u))  =  J  hk(t)eltwdt  =  H*(-uj) 


and  the  noise  spectral  density  matrix  of  the  steered  outputs, 


Gk(")  =  1  RL<'>  dt- 


which  has  the  symmetry  properties 


gl<">  -  Gy-'”>  =  g«^>  • 

In  terms  of  the  spectral  density  matrix,  G  .(uu),  of  the  unsteered  array  outputs,  G*  (uu) 


is  given  by  the  simple  relation 


GMW  ° 


— k  —l' 


(16) 


The  transform  of  the  constraints  (7)  is  simply 


N 

Yj  H  (uj)  =  1  for  all  uu  , 

k=l  * 


(17) 


*  Levin's  h  (i)  equals  our  0  (-i)  and  his  X.  is  our  X(-i). 

X  K  i 
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and  expression  (8)  is  equal  to 


N  00 

o2  =  Z  J  Gu(U))  Hk(UJ)  HI(UU)  duj/2rr 

k,-t=l-“ 


It  can  be  shown  that  G^(ai)  is  a  non-negative  definite  matrix  for  each  UJ,  so  that 
necessarily 


2  Gkf(UU)  Hk(UJ)  “  0  * 


k,«fc=l 


Hence  (18)  is  minimized  by  minimizing  the  integrand,  i.e. ,  the  left  side  of  (19),  at 
each  frequency.  Taking  cognizance  of  the  constraint  (17)  at  each  frequency,  we  obtain 

HkW  =  “n -  (20) 

2  , 

m,n=l 


where  Q'  is  the  inverse  of  the  matrix  G'  at  each  frequency.  The  minimum  noise 
variance  obtained  is  given  by 


o2  =  J*  {  J  Q’  («0 }  du)/2n 


•°°  k,-l=l 


From  (20)  we  find  that  the  H^(uo)  satisfy  the  equations 
N  N 

Z  gm(uj)  Hk(UJ)  =  {  2  q: _ (^)i_1  =  Ar(^) 


m,  n=l 


for  each  l.  In  time  domain,  this  means  that  the  h^(t)  satisfy  the  system  of  integral 
equations 


6 


(23) 


Z  J  Rjc^t-s^  \<s) ds  =  v(t) 

't' 

where 

GO 

X’(t)  =  J  A'(uj)  e  llJJt  duu/2n 

—  co 

We  also  have,  from  (21), 

00 

0s  =  J  A’ (ID)  d(D/2TT  =  \’(0)  (24) 

—  CO 

The  function  -  k'(t)  is  just  the  Lagrangian  multiplier  function  we  mentioned  above  in  a 
direct  time-domain  approach.  Also,  from  (21)  and  (22)  it  is  obvious  that  A’(uu)  is  the 
spectral  density  of  the  noise  component  of  the  output  beam  y(t). 

IQ*  The  Maximum  -Likelihood  Processor 

The  maximum -likelihood  form  of  processing  has  been  discussed  in  detail  in 
Reference  2.  It  is  a  constructive  procedure  for  estimating  both  the  signal  waveform, 
S(t),  and  the  signal  parameters,  a,  considered  to  be  unknown  a  priori.  The  noise  is 
assumed  to  be  gaussian  and  the  estimation  procedure  takes  place  in  two  steps.  First 

/s 

a  tentative  a  is  chosen  and  an  estimate,  S(t;a),  of  the  signal  waveform  is  constructed. 

This  estimate  is,  in  fact,  the  maximum -likelihood  estimate  of  S(t)  for  the  case  in 

which  a  is  known  and  turns  out  to  be  exactly  the  beam,  y’(t),  discussed  in  II,  i.  e. ,  a 

sum  of  the  filtered,  steered  array  outputs  with  filter  functions  given  by  (20). 

In  the  general  case,  a  modified  likelihood  function,  L(cx),  is  also  constructed 

which  is  then  maximized  over  a.  The  resulting  value,  a,  is  the  maximum -likelihood 
*  A 

estimate  of  a  and  S(t;ct)  is  the  final  estimate  of  S(t).  In  terms  of  the  transform,  Y’(uj), 
of  y*(t),  the  function  L(cc)  is  given  by 

dl"/2rr  (24) 
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The  expression,  a  sort  of  signal -to -noise  ratio,  contains  a  implicitly  by  way  of  the 
array  steering  indicated  by  the  primes  in  our  notation.  In  general,  the  final  maximum - 
likelihood  estimate  of  S(t)  is  biased,  as  would  be  any  estimator  for  the  case  of  unknown 
a  and  more  than  one  sensor.  The  accuracy  of  the  estimate  of  a  is  discussed  in  detail 
in  Reference  2. 

It  can  be  shown  that  the  maximum -likelihood  filter  weights  for  the  estimation 
of  signal  with  known  parameters  in  the  sampled -data  case  are  also  identical  with  those 
obtained  by  Levin. 

IV.  Multichannel  Wiener  Filtering 

This  case  was  also  discussed  briefly  in  Reference  2,  but  will  be  reviewed  in 
greater  detail  here.  As  in  II,  we  form  a  steered  beam,  y’(t),  given  by  (6),  but  now 
we  choose  the  filter  functions  to  minimize  the  variance  of  the  error  function 

e(t)  =  y’(t)  -S(t)  ,  (25) 


where  S(t)  is  taken  to  be  a  stationary  random  process  with  covariance  function 


Ro(t-t')  =  E  S(t)  S(t’)  =  Ro(t’-t)  , 

independent  of  the  noise  components. 

First,  we  compute 

E  ^(t)  x^(t*)  =  Ro(t-t')  +  R^(t-t’) 
and 

Ex^(t)S(t’)  =  Ro(t-t’)  , 

hence 
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N 


Ey’(t)2  =  E  SI  hk(T)h^(T’)Ex^(t-T)x^(t-T')drdT’ 

kj  -t — 1 


N 


=  Z  jT  hk(T)h^(T,){R0(T,-T)+R^(T,-T)}  dTdT 


k,^l 


kV 


and 


N 


E  y’(t)  S(t)  =  2  J*  ^(t)  E  x^(t-T)  S(t)  dT 


N 


Z  I  Mt)  r0(_t)  dT 

k=l 


N 


Z  1  Mt)  r0(t)  dT 

k=l 


Finally,  combining  these  results,  we  have 


N 


E{y*(t)-S(t)}  =  2  H  VT)  Vt,)^Ro(t,_t)  +  ru(t,"t)}  dTdT’ 


k,l=l 


N 


"  2  Z  1  hk(T>  R0(t>  ^  +  Rn(0) 


k=l 


At  this  point  we  reintroduce  the  Fourier  transforms  used  in  II  and  also  the 
signal  spectral  density 


Go(W)  =  J*  Ro(t)  elUJt  dt 


Expressing  all  time  functions  in  terms  of  their  Fourier  transforms,  we  find 
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2  N 

E  {y’(t)  -  S(t) }  =  E  J  Hk(iu)  H*(u>)  (Go(uj)  +  G^(oj)}  daj/2n 


N  “ 

~  Z  J  |Hk(UJ)  +  Hj(uJ)}  G  (UI)  duu/2TT 

k=l  —oo 


00 

+  J  Gq(w)  <W2tt  . 
—  00 


We  now  vary  the  Hk(w)  (as  complex  quantities)  and  find  that  the  Wiener  filter  functions 
satisfy  the  equations 


N 

Z  {G^(»)  +  G0(»)} 


Hj>)  =  GqP) 


(26) 


In  terms  of  the  inverse,  QJ^uj),  to  the  matrix  G^(uj),  we  try  a  solution  of  the  form 

N 

Hk(w)  =  H(<D)  Z  Q^k(«J)  Gq(oj)  . 


We  find  that 


H(w) 


{  1  + 


gog)  -1 

A(uj)  1 


and  that  the  Wiener  filter  functions  are 


N 

Z  Q \M 

i»j=l  J 


G0(m) 

G0(ou)  +  A'(uu) 


(27) 


These  functions  are  equal  to  the  MVU  filters,  multiplied  by  the  common  filter  function, 
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G0(u>) 

G0(w)  +  A» 


which  is  simply  the  appropriate  infinite -memory,  non -realizable  Wiener  filter  to  apply 
to  the  MVU  output  beam  to  convert  it  directly  into  a  Wiener  estimate  of  S(t). 

The  smoothing  error  committed  by  the  multichannel  Wiener  filter  is 


N 


E  {y’(t)  —  S(t) }  =  J  Gq(iu)  {  1  —  2  H  (U))}  dUJ/2n 


k=l 


"  G0(ui)  A' (U))  i  _ 
G0(iw)  +  A’ (id)  dU,/2n 


(28) 


The  case  of  unknown  a  is  not  directly  amenable  to  Wiener  smoothing  methods. 

The  sampled  data  version  of  the  Wiener  filtering  problem  is  easily  formulated 

3 

and  formally  solved.  However,  because  of  the  finite  memory  restrictions,  we  shall 
not  find  quite  so  direct  a  relationship  to  the  sampled -data  MVU  solution  as  in  the 
continuous  time  case. 

We  assume  that  the  sampled  signal  sequence,  ^(i)t  is  statistically  stationary 
with  covariance  matrix 

E  £(i)  E(J)  =  P0(i"j)  =  Po(j_i)  ’ 

and  form  a  beam  as  in  IE: 

N  v 

•n(i)  =2  2  Yk(j)  S  (i-j)  •  (29) 

k=l  j=-v 

We  use  Yk(i)  for  the  filter  weights  to  distinguish  them  from  the  MVU  weights  0^.(1).  We 
now  require  a  minimum  value  for 

E  [E(i)-  Z(i)l2  . 

which  leads  immediately  to  the  set  of  equations 


11 


(30) 


N  v 

S  E  [PL.l<j'l)  +  Po(H)1  Yt(J>  =  P0(i) 
tel  j=-v 

for  k  =  1, . . . ,  N  and  i  =  -  v . v. 

In  order  to  solve  this  system  we  modify  the  notation,  and  define  the  doubly- 
indexed  matrices 


*  V<H> 


k,i  =  1 . N 


(31) 


By  A k  we  mean  shall  the  inverse  of  A,  in  the  sense  that 

m=l  s=-v  ’  *  »j 

Equations  (30)  now  read 


N  v 

Z  Z  [Ak  ^i,j) +Bk  =  po(i)  » 

tel  j=-v  ’  ’ 


(32) 


and  our  problem  is  to  invert  the  sum  matrix  A  +  B.  We  wish  to  express  this  inverse 
*  -1  .  -1 

in  terms  of  A  since  A  appears  in  the  MVU  solution.  We  digress  here  to  obtain 
expressions  for  the  MVU  filters  in  terms  of  A"1.  With  our  present  notation,  equation 
(12) reads 

N  v 

E  E  Ak  j)  Mi)  =  -mo  , 

te 1  j=-V 


hence 


N  v 

ek(o  =  -  E  Z  v^ij)  x(j) 

-1=1  j=-v  ' 


(33) 
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According  to  the  constraints,  (10),  we  have 


N  v  N 

6i,o  =  2  9k<‘>  *  -  Z  Z  \  Mi)  • 

k=l  j=-v  k, -t=l 

We  define  the  (2v+l)  x  (2v+l)  matrix  a(i,  j)  by  means  of  its  inverse: 

-v  .  V  -1 

k,tel 


Thus 


V 


6i  o  =  -  2  a"  «.)>  Mi) 
j=-v 


and  we  have 


(34) 


Mi)  ~  ~  2  a(i,j)§  =  -  a(i,o)  (35) 

j=-v  J> 

Finally,  we  have  the  desired  expression  for  the  MVU  filters: 

N  v 

0k^  =  S  S  a(j,o)  (36) 

k=l  j=-v  * 

Had  we  asked  that  the  MVU  beam,  r|(i),  be  an  estimate  of  ^(i-s),  for  some  integer  s 
between  -v  and  +v,  still  using  the  2v+l  data  prints  symmetrically  disposed  about  the 
point  i,  the  analysis  would  be  unchanged,  except  that  in  constraint  (10)  we  have  6.  on 
the  right  side  instead  of  We  call  the  resulting  MVU  filters  0k(i|s),  and  the 

multipliers  that  go  with  them  A(i  |s) ,  and  we  have  immediately 

N  v 

ek(i|s)  =  2  2  Ak  ,(i.j)  a(j,s)  (37) 

tel  j=-v  ’  l'  ’ 

Of  course,  @k(i)  =  0k(i|o)  and  X(i)  =  X(i|o).  We  also  find  that  the  noise  variance  attained 
by  the  filters  ©k(i|s)  is  just  -  X(s  |s)  =  a(s,  s),  while  a(s,  s’)  =  the  covariance  of  the  noise 
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outputs  of  the  two  processors  0^(i|s)  and  0^(i  |s* ).  We  have  introduced  these  modified 
MVU  filters  because  the  Wiener  filter  weights  are  expressed  in  terms  of  a  combination 
of  all  of  them. 

Returning  to  the  Wiener  filters,  we  make  use  of  the  matrix  expansion  analogue 
of 

(l+x)  1  =  Yi  (_1)n  X11  , 
n=o 

namely 

00 

(A+B)"1  =  A_1  J  (-Dn  (BA_1)n  .  (38) 

n=o 

[  Here  (BA  *)°  is  the  unit  matrix .  ] 

Thus 

Yk(l)  =  2  Z  \  ^>3)  Z  (~L)n 

l,m=l  j,s=-v  n=0 


Now  suppose  we  introduce  the  (2v+l)  X  (2v+l)  matrix 


b(i,j)  =  PQ( j-i)  =  Bu(i,j)  . 


Then 


1  j=-v 

‘  Z  f, 

l,  m=l  j,  s=-v 
v  N 

-  Z  b(i,s)  Z  A-*  (s,J)  P0(j) 

j,s=-v  l,  m=l 

v 

=  Z  b(i.s)  a_1(s,j)  pQ(j) 

j»s=-v 

V 

=  Z  [t>a-1]  (i,j)  P0(j)  ,  k=l . N 

i=-v 
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in  terms  of  the  matrix  product  ba-*.  In  fact, 
N  v 

2  2  <i,j)  r0(j) 


£-1  j— “V 


=  J  [ba-l]n  (i,j)P0(j)  ,  k  =  1 . N  , 

j=-v 


and  therefore, 


N  v  “ 

Yk(i)  =  Z  Z  2  (-i)n  [ha_1]n  (j,s)  PQ(j) 

1=1  j=-v  *  n=o 


But 


2  (-l)n  [ba'1]11  =  a(a  +  b)~ 1  , 


n=o 


and  therefore, 

N  v 

Yk(0  =  Z  Z  Ak  a(j,s)  (a+b)_1(s,t)  PQ(j) 
^lj,s,t=-v  ’ 

We  define 


(39) 


W(i)  =  Z  (a  +  b)~ 1  (i,j)  PQ(j) 
j=-v 

and  find  the  desired  expression  in  terms  of  the  MVU  filters: 

v 

Yk(i)  =  S  w(j)  0k(i|j)  (40) 

j=-V 

which  is  analogous  to  equation  (27)  in  the  continuous  case.  It  might  be  remarked  that 
the  sequence  W(i)  is  vaguely  related  to  the  single -channel  "signal /signal  noise"  Wiener 
filter.  Finally,  the  Wiener  smoothing  error  is 
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(41) 


E  [n(i)  -  2(i)]2  =  Po(0)  -  t  (a  +  b)" 1  (i,j)  PQ(i)  P0(j) 

We  can  also  introduce  a  set  of  Wiener  filter  weights,  y^(i|s),  in  analogy  to  the 
MVU  filters  9^(i|s),  by  defining 

N  v 

Is)  =  E  E  Yv(ils)  ?’(i-j)  (42) 

k=l  j=-v  k  k 

to  be  a  Wiener  estimate  of  E(i_s)-  The  only  change  is  that  PQ(i)  becomes  PQ(i-s)  on 
the  right  side  of  (30)  and  W(j)  in  (40)  becomes  W(j|s),  where 


v 


W(i|s)  =  E  (a+b)  (i,j)P0(i-s) 

j=-v 


(43) 


The  corresponding  smoothing  error  is 


E  [T!(i|s)  -  2(i-s)]2  =  Po(0)“  E  (a  +  b)_1(i,j)P0(i-s)Po(j-s) 

i»j=-v 

=  {b  —  b(a  +  b)"-*-  b}  (s,s) 

=  {a  (a  +b)  1  b}  (s,s)  (44) 


Incidentally, 

v 

W(i|s)  =  E  (a  +  b)  1  (i,  j)  b(j,s) 
j=-V 

=  {  (a  +b)  1  b}  (i,  s)  . 
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